例如一位工程师想找出对桥梁恶化最重要对影响因素(时长,桥梁设计,材料,天气等),来确定这些关系的数学形式。以下几点会有助于这个过程:
本章的思路: 1. 如何表示模型 2. 展示各种回归方法 3. 判断回归模型是不是有问题 4. 找出特别的观测值 5. 纠正模型
如何用 R 表示 输入与响应的公式? Y ~ X1 + X2 + ... + Xk 就表示一个线性回归的公式。
| 符号 | 应用 |
|---|---|
~ |
用来分开输入与响应,左边的位响应 |
+ |
间隔开输入变量 |
: |
两个输入变量间的相互作用 , y ~ x + z + x:z |
* |
所有有可能的相互作用的简写, y ~ x * z * w <=> y ~ x + z + w + x:z + x:w + z:w + x:z:w |
^ |
标记相互作用的上限, y ~ (x + z + w)^2 <=> y ~ x + z + w + x:z + x:w + z:w |
. |
data frame 中除去依赖变量的所有变量的占位符, y ~ . <=> y ~ x + z + w |
- |
从方程中移除某项, y ~ (x + z + w)^2 - x:w <=> y ~ x + z + w + x:z + z:w |
-1 |
去掉截距,例如: y ~ x -1 强制回归线通过原点 |
I() |
表示括号中的需要预先处理,例如:y ~ x + I((z + w)^2) 将会扩展成 y ~ x + h,h 是一个新变量,\(h=\sum(z+w)^2\) |
function |
将数学函数使用在方程中, log(y) ~ x + z + w,将预测 log(y) 而不是 y |
\[ RSS=\sum^n_{i=1}(y_i - \hat{y_i})^2 \]
fit <- lm(weight ~ height, data=women)
plot(women$height, women$weight)
abline(fit)
summary(fit)
##
## Call:
## lm(formula = weight ~ height, data = women)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7333 -1.1333 -0.3833 0.7417 3.1167
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -87.51667 5.93694 -14.74 1.71e-09 ***
## height 3.45000 0.09114 37.85 1.09e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.525 on 13 degrees of freedom
## Multiple R-squared: 0.991, Adjusted R-squared: 0.9903
## F-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14
residuals(fit)
## 1 2 3 4 5 6
## 2.41666667 0.96666667 0.51666667 0.06666667 -0.38333333 -0.83333333
## 7 8 9 10 11 12
## -1.28333333 -1.73333333 -1.18333333 -1.63333333 -1.08333333 -0.53333333
## 13 14 15
## 0.01666667 1.56666667 3.11666667
fitted(fit)
## 1 2 3 4 5 6 7 8
## 112.5833 116.0333 119.4833 122.9333 126.3833 129.8333 133.2833 136.7333
## 9 10 11 12 13 14 15
## 140.1833 143.6333 147.0833 150.5333 153.9833 157.4333 160.8833
plot(fit)
anova(fit)
## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## height 1 3332.7 3332.7 1433 1.091e-14 ***
## Residuals 13 30.2 2.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
vcov(fit)
## (Intercept) height
## (Intercept) 35.247305 -0.539880952
## height -0.539881 0.008305861
回归系数为3.45显著不同于0(p < 0.001)
fit2 <- lm(weight ~ height + I(height^2), data=women)
summary(fit2)
##
## Call:
## lm(formula = weight ~ height + I(height^2), data = women)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50941 -0.29611 -0.00941 0.28615 0.59706
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 261.87818 25.19677 10.393 2.36e-07 ***
## height -7.34832 0.77769 -9.449 6.58e-07 ***
## I(height^2) 0.08306 0.00598 13.891 9.32e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3841 on 12 degrees of freedom
## Multiple R-squared: 0.9995, Adjusted R-squared: 0.9994
## F-statistic: 1.139e+04 on 2 and 12 DF, p-value: < 2.2e-16
plot(women$height, women$weight)
lines(women$height,fitted(fit2))
library(car)
scatterplot(weight ~ height, data=women, spread=FALSE, lty.smooth=2, pch=19, main="Women Age 30-39", xlab="Height (inches)", ylab="Weight (lbs.)")
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
cor(states)
## Murder Population Illiteracy Income Frost
## Murder 1.0000000 0.3436428 0.7029752 -0.2300776 -0.5388834
## Population 0.3436428 1.0000000 0.1076224 0.2082276 -0.3321525
## Illiteracy 0.7029752 0.1076224 1.0000000 -0.4370752 -0.6719470
## Income -0.2300776 0.2082276 -0.4370752 1.0000000 0.2262822
## Frost -0.5388834 -0.3321525 -0.6719470 0.2262822 1.0000000
scatterplotMatrix(states, spread=FALSE, main="Scatter Plot Matrix")
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
summary(fit)
##
## Call:
## lm(formula = Murder ~ Population + Illiteracy + Income + Frost,
## data = states)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7960 -1.6495 -0.0811 1.4815 7.6210
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.235e+00 3.866e+00 0.319 0.7510
## Population 2.237e-04 9.052e-05 2.471 0.0173 *
## Illiteracy 4.143e+00 8.744e-01 4.738 2.19e-05 ***
## Income 6.442e-05 6.837e-04 0.094 0.9253
## Frost 5.813e-04 1.005e-02 0.058 0.9541
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.535 on 45 degrees of freedom
## Multiple R-squared: 0.567, Adjusted R-squared: 0.5285
## F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08
library(effects)
fit <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)
summary(fit)
##
## Call:
## lm(formula = mpg ~ hp + wt + hp:wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0632 -1.6491 -0.7362 1.4211 4.5513
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.80842 3.60516 13.816 5.01e-14 ***
## hp -0.12010 0.02470 -4.863 4.04e-05 ***
## wt -8.21662 1.26971 -6.471 5.20e-07 ***
## hp:wt 0.02785 0.00742 3.753 0.000811 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.153 on 28 degrees of freedom
## Multiple R-squared: 0.8848, Adjusted R-squared: 0.8724
## F-statistic: 71.66 on 3 and 28 DF, p-value: 2.981e-13
plot(effect("hp:wt", fit, vcov, list(wt=c(2.2,3.2,4.2))), multiline=TRUE)
上述的 lm 方法,只会根据我们给的公式进行拟合,并产生模型,并没有输出说模型是不是合适。有时还会发生,明明没有关系的输入与响应,最后的结论却是二者有关。
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
confint(fit)
## 2.5 % 97.5 %
## (Intercept) -6.552191e+00 9.0213182149
## Population 4.136397e-05 0.0004059867
## Illiteracy 2.381799e+00 5.9038743192
## Income -1.312611e-03 0.0014414600
## Frost -1.966781e-02 0.0208304170
以上结果表明,文盲率1%的改变,谋杀率有95%的可能在 \([2.38, 5.90]\) 间。因为 Frost 置信区间包含0,我们可以得出谋杀率和气温没有关系。这些推断只是在数据符合模型下的统计假设时才的出的。
fit <- lm(weight ~ height, data=women)
par(mfrow=c(2,2))
plot(fit)
fit2 <- lm(weight ~ height + I(height^2), data=women)
par(mfrow=c(2,2))
plot(fit2)
# 去掉异常值
newfit <- lm(weight~ height + I(height^2), data=women[-c(13,15),])
par(mfrow=c(2,2))
plot(newfit)
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
par(mfrow=c(2,2))
plot(fit)
library(car)
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
qqPlot(fit, labels=row.names(states), simulate=TRUE, main="Q-Q Plot")
states["Nevada",]
## Murder Population Illiteracy Income Frost
## Nevada 11.5 590 0.5 5149 188
fitted(fit)["Nevada"]
## Nevada
## 3.878958
residuals(fit)["Nevada"]
## Nevada
## 7.621042
rstudent(fit)["Nevada"]
## Nevada
## 3.542929
从上我们看出,Nevada 有很高的谋杀率(11.5),但是预测出的值确很小(3.878958)。问题出在哪里?
residplot <- function(fit, nbreaks = 10) {
z <- rstudent(fit)
hist(z, breaks = nbreaks, freq = FALSE, xlab="Studentized Residual",
main="Distribution of Errors")
rug(jitter(z), col = "brown")
curve(dnorm(x, mean = mean(z), sd = sd(z)), add = TRUE,
col = "blue", lwd = 2)
lines(density(z)$x, density(z)$y, col = "red", lwd = 2, lty = 2)
legend("topright",
legend = c( "Normal Curve", "Kernel Density Curve"),
lty=1:2, col=c("blue","red"), cex=.7)
}
residplot(fit)
rstudent(fit)["Nevada"]
## Nevada
## 3.542929
误差应当符合正态分布。
可以注意到右下角出现了异常值。
# 判断是否有自相关性
durbinWatsonTest(fit)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.2006929 2.317691 0.248
## Alternative hypothesis: rho != 0
# lag为1,表明每个观测值都是和下一个值进行比较的。
# 判断是否有线性关系
crPlots(fit)
# 同方差性
ncvTest(fit)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 1.746514 Df = 1 p = 0.1863156
spreadLevelPlot(fit)
##
## Suggested power transformation: 1.209626
推荐的指数 p(\(Y^p\))
library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)
##
## Call:
## lm(formula = Murder ~ Population + Illiteracy + Income + Frost,
## data = states)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7960 -1.6495 -0.0811 1.4815 7.6210
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.235e+00 3.866e+00 0.319 0.7510
## Population 2.237e-04 9.052e-05 2.471 0.0173 *
## Illiteracy 4.143e+00 8.744e-01 4.738 2.19e-05 ***
## Income 6.442e-05 6.837e-04 0.094 0.9253
## Frost 5.813e-04 1.005e-02 0.058 0.9541
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.535 on 45 degrees of freedom
## Multiple R-squared: 0.567, Adjusted R-squared: 0.5285
## F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = fit)
##
## Value p-value Decision
## Global Stat 2.7728 0.5965 Assumptions acceptable.
## Skewness 1.5374 0.2150 Assumptions acceptable.
## Kurtosis 0.6376 0.4246 Assumptions acceptable.
## Link Function 0.1154 0.7341 Assumptions acceptable.
## Heteroscedasticity 0.4824 0.4873 Assumptions acceptable.
多重共线性是:多个输入变量组合起来的时候和响应有关系,但是单个输入变量确和响应没有关系。
vif(fit)
## Population Illiteracy Income Frost
## 1.245282 2.165848 1.345822 2.082547
sqrt(vif(fit)) > 2 # problem?
## Population Illiteracy Income Frost
## FALSE FALSE FALSE FALSE
以上检测说明我们的输入变量中不存在多重共线性的问题。
异常值是指那些预测值与实际值偏差很大的值。
outlierTest(fit)
## rstudent unadjusted p-value Bonferonni p
## Nevada 3.542929 0.00095088 0.047544
高杠杆率点:
high leverage are outliers with regard to the other predictors. In other words, they have an unusual combination of predictor values. The response value isn’t involved in determining leverage.
高杠杆率点可以通过 hat 统计的方法来得到,对于给定的数据集,hat 值的平均值为 \(p/n\) ,其中 p 为模型中估计参数的个数,包括截距,n 为数据集的条目数。大于这个平均值 2 到 3 倍的观测值的 hat 值,都是高杠杆率点。
hat.plot <- function(fit) {
p <- length(coefficients(fit))
n <- length(fitted(fit))
plot(hatvalues(fit), main="Index Plot of Hat Values")
abline(h = c(2,3)*p/n, col = "red", lty = 2)
identify(1:n, hatvalues(fit), names(hatvalues(fit)))
}
hat.plot(fit)
## integer(0)
影响点是那些对模型参数有不成比例的影响的观测值。例如,去掉某个观测值,模型会发生很大的变动。
Cook’s D 值大于 \(4/(n-k-1)\), k 是输入变量的个数, 就表明是影响点。
cutoff <- 4/(nrow(states)-length(fit$coefficients)-2)
plot(fit, which=4, cook.levels=cutoff)
abline(h=cutoff, lty=2, col="red")
有时候 Cook’s D 为1,比 \(4/(n-k-1)\) 更有用。
Cook’s D 只能标识出影响点,却没有这些点怎样影响模型的信息。 Added-Variable Plots 可以给出这些信息。
avPlots(fit, ask=FALSE, id.method="identify")
coefficients(fit)
## (Intercept) Population Illiteracy Income Frost
## 1.2345634112 0.0002236754 4.1428365903 0.0000644247 0.0005813055
每幅图中的直线表示该输入变量的系数。
influencePlot 可以画出以上提到的3中点
influencePlot(fit, id.method="identify", main="Influence Plot",
sub="Circle size is proportional to Cook’s distance")
在水平线 +2,-2 上下的点是异常值,垂直线0.2或 0.3右边的点。圆圈的大小表示该观测值对模型影响的大小。
通常是删掉影响点或者异常值,然后再生成模型。这种删除要慎重。
当模型违反正态假设的时候,一种是修改响应变量的次数,\(Y\rightarrow Y^\lambda\)。如果 \(Y\) 是一个比值,则可以用 \(ln(Y/1-Y)\)
summary(powerTransform(states$Murder))
## bcPower Transformation to Normality
##
## Est.Power Std.Err. Wald Lower Bound Wald Upper Bound
## states$Murder 0.6055 0.2639 0.0884 1.1227
##
## Likelihood ratio tests about transformation parameters
## LRT df pval
## LR test, lambda = (0) 5.665991 1 0.01729694
## LR test, lambda = (1) 2.122763 1 0.14512456
从上可以看出,推荐的参数为0.6(\(Murder^{0.6}\)),因为0.6接近与0.5,则可以考虑使用开平方。但是在这个例子中 \(\lambda=1\) 不能被拒绝掉(\(p=0.145\)),则没有强有力的证据,证明需要这样一个变换。
另一种是对输入变量进行变换,
boxTidwell(Murder~Population+Illiteracy,data=states)
## Score Statistic p-value MLE of lambda
## Population -0.3228003 0.7468465 0.8693882
## Illiteracy 0.6193814 0.5356651 1.3581188
##
## iterations = 19
以上结果表明 \(Population^{0.87}\,,Population^{1.36}\)可以提升线性关系。但二者的 p 值表明没必要进行这种变换。
要处理多重共线性,删除变量是一个很重要的方法。当我们只是进行预测,多重共线性对我们没有什么影响。而当我们要理解预测值与每个输入之间对关系时,就需要处理多重共线性。一种方法时删除 \(\sqrt{vif}>2\) 对变量。另一种是使用脊回归。
选择精准的模型,选择简单的模型。
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
fit2 <- lm(Murder ~ Population + Illiteracy, data=states)
anova(fit1)
## Analysis of Variance Table
##
## Response: Murder
## Df Sum Sq Mean Sq F value Pr(>F)
## Population 1 78.854 78.854 12.2713 0.001052 **
## Illiteracy 1 299.646 299.646 46.6307 1.83e-08 ***
## Income 1 0.057 0.057 0.0089 0.925368
## Frost 1 0.021 0.021 0.0033 0.954148
## Residuals 45 289.167 6.426
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fit2)
## Analysis of Variance Table
##
## Response: Murder
## Df Sum Sq Mean Sq F value Pr(>F)
## Population 1 78.854 78.854 12.813 0.0008122 ***
## Illiteracy 1 299.646 299.646 48.690 8.832e-09 ***
## Residuals 47 289.246 6.154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fit2, fit1)
## Analysis of Variance Table
##
## Model 1: Murder ~ Population + Illiteracy
## Model 2: Murder ~ Population + Illiteracy + Income + Frost
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 47 289.25
## 2 45 289.17 2 0.078505 0.0061 0.9939
模型1包含在模型2中,因为非显著性测试为0.9939,则说明二者相同,所以我们第2个模型,即去掉 Income 和 Frost。
Akaike Information Criterion (AIC) 提供了另一种比较模型的方法(考虑进了达到拟合的参数个数)。更小的 AIC 值,意味着较少的参数就达到了较高的精准度,也更适合。
AIC(fit1,fit2)
## df AIC
## fit1 6 241.6429
## fit2 4 237.6565
需要注意的是 anova 需要嵌套的模型,而 AIC 则不需要。
添加/删除变量,直到达到某个标准。例如前向逐步回归是每次增加一个变量,直到模型不再提升;后向逐步回归则是一开是包含所有变量,每次去掉一个变量,直到模型开始降低。逐步回归则是前向+后向。
library(MASS)
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
stepAIC(fit, direction="backward")
## Start: AIC=97.75
## Murder ~ Population + Illiteracy + Income + Frost
##
## Df Sum of Sq RSS AIC
## - Frost 1 0.021 289.19 95.753
## - Income 1 0.057 289.22 95.759
## <none> 289.17 97.749
## - Population 1 39.238 328.41 102.111
## - Illiteracy 1 144.264 433.43 115.986
##
## Step: AIC=95.75
## Murder ~ Population + Illiteracy + Income
##
## Df Sum of Sq RSS AIC
## - Income 1 0.057 289.25 93.763
## <none> 289.19 95.753
## - Population 1 43.658 332.85 100.783
## - Illiteracy 1 236.196 525.38 123.605
##
## Step: AIC=93.76
## Murder ~ Population + Illiteracy
##
## Df Sum of Sq RSS AIC
## <none> 289.25 93.763
## - Population 1 48.517 337.76 99.516
## - Illiteracy 1 299.646 588.89 127.311
##
## Call:
## lm(formula = Murder ~ Population + Illiteracy, data = states)
##
## Coefficients:
## (Intercept) Population Illiteracy
## 1.6515497 0.0002242 4.0807366
<none> 的 AIC 值表示没有变量移除时,模型的 AIC。
第一步中,Frost 被移除,AIC 从97.75下降到95.75。该方法有个限制就是并不会评估所有的模型。而下面要介绍的这个方法讲会克服这个限制。
library(leaps)
leaps <-regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data=states, nbest=4)
plot(leaps, scale="adjr2")
上图显示,最好的模型是两个变量的模型(Population 和 Illiteracy)。
subsets(leaps, statistic="cp", main="Cp Plot for All Subsets Regression", legend=c(3.5, 50))
abline(1,1,lty=2,col="red")
k-fold 交叉验证:将样本分成 k 个子样本,每个子样本作为测试集,同时除了本身之外的所有子集会构成训练集。
library(bootstrap)
shrinkage <- function(fit, k=10){
theta.fit <- function(x,y){ lsfit(x,y) }
theta.predict <- function(fit,x){ cbind(1,x) %*% fit$coef }
x <- fit$model[,2:ncol(fit$model)]
y <- fit$model[,1]
results <- crossval(x, y, theta.fit, theta.predict, ngroup=k)
r2 <- cor(y, fit$fitted.values)^2
r2cv <- cor(y, results$cv.fit)^2
cat("Original R-square =", r2, "\n")
cat(k, "Fold Cross-Validated R-square =", r2cv, "\n")
cat("Change =", r2-r2cv, "\n")
}
fit <- lm(Murder ~ Population + Income + Illiteracy + Frost, data=states)
shrinkage(fit)
## Original R-square = 0.5669502
## 10 Fold Cross-Validated R-square = 0.494881
## Change = 0.07206925
fit2 <- lm(Murder~Population+Illiteracy,data=states)
shrinkage(fit2)
## Original R-square = 0.5668327
## 10 Fold Cross-Validated R-square = 0.5281635
## Change = 0.03866918
哪个变量在预测响应的时候最重要?
一种是看标准回归系数的大小
zstates <- as.data.frame(scale(states))
zfit <- lm(Murder~Population + Income + Illiteracy + Frost, data=zstates)
coef(zfit)
## (Intercept) Population Income Illiteracy Frost
## -2.054026e-16 2.705095e-01 1.072372e-02 6.840496e-01 8.185407e-03
再有就是 relaimpo 包中的方法。
relweights <- function(fit,...){
R <- cor(fit$model)
nvar <- ncol(R)
rxx <- R[2:nvar, 2:nvar]
rxy <- R[2:nvar, 1]
svd <- eigen(rxx)
evec <- svd$vectors
ev <- svd$values
delta <- diag(sqrt(ev))
lambda <- evec %*% delta %*% t(evec)
lambdasq <- lambda ^ 2
beta <- solve(lambda) %*% rxy
rsquare <- colSums(beta ^ 2)
rawwgt <- lambdasq %*% beta ^ 2
import <- (rawwgt / rsquare) * 100
lbls <- names(fit$model[2:nvar])
rownames(import) <- lbls
colnames(import) <- "Weights"
barplot(t(import),names.arg=lbls,
ylab="% of R-Square",
xlab="Predictor Variables",
main="Relative Importance of Predictor Variables",
sub=paste("R-Square=", round(rsquare, digits=3)),
...)
return(import)
}
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
relweights(fit, col="lightgrey")
## Weights
## Population 14.723401
## Illiteracy 59.000195
## Income 5.488962
## Frost 20.787442